NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIFORNIA 


THESIS 


AUTOMATING  IDENTIFICATION  OF  ROADS  AND 
TRAILS  UNDER  CANOPY  USING  LIDAR 

by 

Charles  F.  Harmon  III 
September  2011 

Thesis  Advisor:  Richard  C.  Olsen 

Second  Reader:  Kristen  Tsolis 


Approved  for  public  release;  distribution  is  unlimited 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to 
Washington  headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA 
22202^1302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank )  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

September  2011  Master’s  Thesis 


4.  TITLE  AND  SUBTITLE  Automating  Identification  of  Roads  and  Trails  Under  5.  FUNDING  NUMBERS 
Canopy  Using  LID  AR 
6.  AUTHOR  Charles  F.  Harmon  III 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 


9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy 
or  position  of  the  Department  of  Defense  or  the  U.S.  Government.  IRB  Protocol  number  NA. 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT  12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  is  unlimited 


13.  ABSTRACT  (maximum  200  words) 

Analysis  techniques  are  developed  to  automatically  extract  roads  and  trails  under  thick  forest  canopy.  LiDAR  data 
were  taken  over  the  Swanton  Pacific  Ranch  in  the  Santa  Cruz  Mountains  from  an  airborne  laser  mapping  system,  the 
Optech  3100,  on  March  9—10,  2010.  Collected  data  were  characterized  by  point  densities  of  5—10  m2.  Point  cloud 
data  were  reduced  to  digital  surface  models  using  ARCMAP  (from  ESRI).  The  DSM  was  calculated  at  1  meter 
spacing.  These  surface  models  were  analyzed  using  topographic  tools  in  EN  VI,  allowing  for  calculation  of  curvature, 
slope,  convexity,  and  shaded  relief.  A  multi-layer  dataset  was  built  and  analyzed  using  spectral  analysis  tools  in 
ENVI.  The  classification  technique  used  was  a  combination  of  maximum  likelihood  classifier  and  a  decision  tree 
after  use  of  erosion/dilation  operators.  Results  are  compared  to  ground  truth  collected  in  2011.  Classification 
resulted  in  83.6%  true  positive  rate,  and  the  image  processing  result  reduced  the  false  positive  rate  to  3.0%. 


14.  SUBJECT  TERMS  LiDAR,  LADAR,  Swanton  Pacific  Ranch,  Feature  Extraction 


17.  SECURITY 
CLASSIFICATION  OF 
REPORT 

Unclassified 


18.  SECURITY 
CLASSIFICATION  OF  THIS 
PAGE 

Unclassified 


19.  SECURITY 
CLASSIFICATION  OF 
ABSTRACT 

Unclassified 


15.  NUMBER  OF 
PAGES 

91 


16.  PRICE  CODE 


20.  LIMITATION  OF 
ABSTRACT 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release,  distribution  is  unlimited 


AUTOMATING  IDENTIFICATION  OF  ROADS  AND  TRAILS  UNDER 

CANOPY  USING  LIDAR 


Charles  F.  Hannon  III 
Major,  United  States  Army 

B.A.,  University  of  North  Carolina  -  Chapel  Hill,  1999 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  SPACE  SYSTEMS  OPERATIONS 

and 

MASTER  OF  SCIENCE  IN  REMOTE  SENSING  INTELLIGENCE 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
September  2011 

Author:  Charles  F.  Harmon  III 


Approved  by:  Richard  C.  Olsen 

Thesis  Advisor 


Kristen  Tsolis 
Second  Reader 


Dan  C.  Boger 

Chair,  Department  of  Information  Science 


Rudy  Panholzer 

Chair,  Department  of  Space  Systems  Academic  Group 
iii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 
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Mountains  from  an  airborne  laser  mapping  system,  the  Optech  3100,  on  March  9-10, 
2010.  Collected  data  were  characterized  by  point  densities  of  5-10  m  .  Point  cloud 
data  were  reduced  to  digital  surface  models  using  ARCMAP  (from  ESRI).  The  DSM 
was  calculated  at  1  meter  spacing.  These  surface  models  were  analyzed  using 
topographic  tools  in  ENVI,  allowing  for  calculation  of  curvature,  slope,  convexity,  and 
shaded  relief.  A  multi-layer  dataset  was  built  and  analyzed  using  spectral  analysis  tools 
in  ENVI.  The  classification  technique  used  was  a  combination  of  maximum  likelihood 
classifier  and  a  decision  tree  after  use  of  erosion/dilation  operators.  Results  are 
compared  to  ground  truth  collected  in  2011.  Classification  resulted  in  83.6%  true 
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I.  INTRODUCTION 


A.  PURPOSE  OF  RESEARCH 

During  intelligence  preparation  of  the  battlefield  (IPB),  intelligence  analysts  must 
describe  the  operational  environment  for  the  commander.  This  includes  mapping  the 
lines  of  communications  (LOC)  that  traverse  the  area  of  operations.  Traditional  methods 
of  remote  sensing  cannot  identify  roads  and  trails  beneath  dense  canopy.  LiDAR  (Light 
Detection  and  Ranging)  is  a  remote  sensing  technology  used  to  create  digital  elevation 
models  (DEM)  of  the  Earth’s  surface,  and  has  the  capability  to  penetrate  forest  canopy  to 
identify  roads  and  trails  beneath.  This  research  will  use  classification  methods  and  image 
processing  techniques  to  automatically  identify  trails  beneath  dense  forest  canopy,  to 
provide  a  tool  for  intelligence  analysts  in  support  of  ground  maneuver  forces. 

LiDAR,  like  RADAR,  is  an  acronym  which  stands  for  light  detection  and  ranging, 
and  also  describes  the  process  by  which  LiDAR  systems  create  3D  point  cloud  models  of 
terrain  and  above  ground  objects  (Gordon  &  Charles,  2008).  LiDAR  systems  measure 
the  time-of-flight  for  a  laser  pulse  to  travel  from  a  sensor  to  a  reflective  object  and  back. 
Over  the  last  decade,  topographic  laser  profiling  and  scanning  (LiDAR)  systems  have 
made  major  improvements  in  accuracy  and  application.  Installed  on  airborne, 
spacebome,  or  terrestrial-based  platforms  these  sensors  are  approaching  horizontal  and 
vertical  accuracies  on  the  order  of  centimeters  (Vossehnan  &  Maas,  2010).  These 
sensors  fire  a  beam  of  light  that  travels  toward  the  target,  i.e.,  ground.  If  the  beam  strikes 
anything  on  its  way  to  the  ground,  such  as  a  tree,  part  of  the  beam’s  energy  is  reflected 
back  to  the  sensor  and  recorded  as  the  first  return.  The  rest  of  the  beam  continues  to  the 
ground  or  some  other  solid  surface  that  prevents  further  progress  and  reflects  back  to  the 
sensor  as  the  last  return  (Crutchley  &  Crow,  2009).  Early  systems  generally  only 
recorded  the  first  and  last  returns  since  these  measurements  were  used  to  create  the  digital 
surface  models  (DSM)  and  digital  terrain  model  (DTM)  respectively.  Modem  systems 
record  up  to  an  industry  standard  four  returns,  and  a  few  systems  are  capable  of  providing 
full  wavefonn  data  as  described  in  Figure  1 . 
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Beam 

divergence 


Figure  1.  Full  waveform  LiDAR  versus  discrete  recording  characteristics  (From  Diaz, 
2011) 


Of  the  total  land  mass  of  Earth,  forests  cover  31%;  this  was  estimated  to  be  over 
10  billion  acres  in  2010  (Food  and  Agriculture  Organization  of  the  United  Nations, 
2010).  Traditional  imaging  systems  are  unable  to  observe  beneath  these  forest  canopies, 
allowing  freedom  of  maneuver  for  illicit  organizations  and  terrorist  activities.  As 
discussed  earlier,  LiDAR  can  penetrate  forest  canopies  and  identify  the  surface 
characteristics  and  any  manmade  structures  under  cover  of  trees.  Previous  work  on  this 
subject  has  proven  the  feasibility  of  using  LiDAR  to  identify  roads  and  trails  beneath 


dense  canopy  (Espinoza  &  Owens,  2007).  The  purpose  of  this  research  is  to  develop  an 
automated  process  to  identify  roads  and  trails  under  tree  cover  from  large  LiDAR  data 
sets. 

B.  OBJECTIVE 

The  primary  objective  of  this  thesis  is  to  detect  roads  and  trails  in  a  forested 
region  of  Santa  Cruz  County,  CA,  from  airborne  LiDAR  data.  The  LiDAR  data  set  was 
collected  in  2010  for  California  Polytechnic  State  University  (Cal  Poly)  for  their  Swanton 
Pacific  Ranch  study  area.  Using  products  derived  from  the  LiDAR  point  cloud,  roads 
and  trails  are  characterized  and  classified  using  a  decision  tree  approach,  where  the  union 
of  a  maximum  likelihood  classification  and  a  Laplacian  edge  enhancement  produces  the 
final  product. 

A  brief  history  of  the  development  of  LiDAR  is  presented  in  the  background 
section  as  well  as  a  discussion  of  the  different  terrestrial,  airborne,  and  spaceborne 
applications  of  LiDAR.  A  short  explanation  of  laser  fundamentals  and  laser  ranging  is 
included  to  illustrate  the  basic  measurement  principles  of  LiDAR  systems.  The  problem 
and  observations  sections  describe  the  areas  of  study,  the  equipment  and  software  used, 
and  the  experimental  setup.  The  analysis  and  summary  sections  provide  conclusions 
drawn  from  the  research  and  experimental  results. 
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II.  BACKGROUND 


A.  BACKGROUND 

1.  Brief  History  of  LiDAR  Development 

RADAR  (Radio  Detection  And  Ranging)  is  the  process  of  transmitting,  receiving, 
detecting,  and  processing  electromagnetic  waves,  and  was  invented  by  the  Gennan  Army 
in  1935  (Richmond  &  Cain,  2010).  Since  then,  this  process  has  expanded  to  use  in  other 
areas  of  the  magnetic  spectrum  from  centimeter  waves,  to  millimeter  and  microwave 
ranges,  and  with  the  advent  of  lasers  to  optical  wavelengths. 

In  1960,  the  invention  of  the  ruby  laser  was  followed  by  a  decade  of  rapid 
development  for  laser  technology.  Shortly  after,  laser  surveying  instruments  were  proven 
in  experimental  use  with  the  first  laser  altimeters  flown  from  aircraft  as  early  as  1965  in 
the  United  Kingdom  by  Shepherd.  Later  that  same  year,  Miller,  Jensen  and  Ruddock 
introduced  the  first  airborne  laser  profiler  for  commercial  topographic  mapping  by  a  joint 
venture  of  Spectra  Physics  Company  (built  the  laser)  and  Aero  Service  Corporation 
(aerial  survey  and  mapping  company)  (Gordon  &  Charles,  2008).  During  the  1970s  and 
1980s,  laser  profiling  systems  experienced  steady  development,  typically  using  Nd:YAG 
solid  state  lasers  or  GaAs  semiconductor  lasers.  The  systems  developed  during  this  time 
used  two-axis  gyros  to  measure  the  aircraft  attitude  and  used  microwave  transponders 
that  calculated  aircraft  position  and  altitude  based  on  triangulation  from  three  surveyed 
ground  stations. 

It  was  not  until  the  advent  of  global  positioning  system  (GPS),  inertial  measuring 
units  (IMU)  for  aircraft,  and  improved  computer  processing  in  the  1990s  that  laser 
scanners  could  achieve  the  accuracies  needed  to  make  them  commercially  viable  for 
topographic  mapping  (Heritage  &  Large,  2009).  In  1989,  Dr.  Joachim  Lindenberger 
conducted  extensive  testing  of  an  airborne  laser  profiling  system  (ALPS)  that  used  a 
Sercel  GPS  receiver  in  conjunction  with  a  Delco  Carrousel  IMU  for  position  and  attitude 
detennination.  Optech  added  a  scanning  mechanism  to  the  system  three  years  later.  The 
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resulting  Airborne  Laser  Terrain  Mapper  (ALTM)  1020  was  the  first  LiDAR  system  that 
had  all  of  the  components  of  modem  systems  (Gordon  &  Charles,  2008). 


Figure  2.  Basic  operating  principles  of  airborne  mapping  LiDAR  with  enabling 

technologies:  Global  Navigation  Satellite  System  (GNSS),  IMU,  and  laser 
scanner  (From  Diaz,  2011) 

2.  Laser  Fundamentals 

Within  the  context  of  LiDAR,  any  point  on  the  Earth’s  surface  can  be  described 
by  its  x,  y,  and  z  coordinates.  A  LiDAR  system  can  construct  these  points  from  three 
sources:  the  LiDAR  sensor,  the  IMU  of  the  aircraft,  and  the  GPS  (Heritage  &  Large, 
2009).  A  scanning  LiDAR  system  corrects  measurements  for  the  yaw,  pitch,  and  roll  of 
the  aircraft  and  the  side-to-side  scanning  mechanism  of  the  LiDAR  sensor.  The  GPS 
units  allow  the  measurements  to  be  overlaid  relative  to  the  WGS84  datum.  The 
following  discussion  describes  lasers  in  general  and  laser  ranging  in  detail. 

LiDAR  shares  many  of  the  same  characteristics  as  radar,  such  as  waveform  and 
propagation  time,  albeit  at  a  different  frequency  band.  Where  radars  operate  over  a  wide 
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range  of  frequencies,  usually  based  on  target  characteristics  and  atmospheric  attenuation, 
lasers  are  limited  to  discrete  laser  lines  and  usually  operate  in  the  near  infrared  based  on 
acceptable  lasers  and  detectors  (Richmond  &  Cain,  2010).  LiDAR  systems’  pulse 
repetition  frequency  (PRF)  greatly  affects  resolution;  as  systems  with  higher  PRFs  can 
measure  more  points  in  a  given  area  of  study. 

Generally,  lasers  are  classified  by  the  material  used  as  a  radiation  source,  most 
commonly  gas  lasers,  solid-state  lasers,  and  semi-conductor  lasers.  Lasers  used  for 
topographic  mapping  are  required  to  have  high  intensity  and  have  a  high  degree  of 
collimation,  that  is  the  light  rays  are  near  parallel,  and  will  spread  slowly  as  it  propagates 
(Gordon  &  Charles,  2008).  Given  these  requirements,  the  most  common  laser  types  used 
for  topographic  mapping  are  solid-state  lasers  using  neodymium-doped  yttrium 
aluminum  garnet  (Nd:YAG)  and  semiconductor  lasers  using  gallium  arsenide  (GaAs). 
These  materials  when  coupled  with  an  energy  source  and  two  mirrors,  one  fully 
reflective,  and  the  other  semi  reflective  make  up  the  components  of  every  laser. 


Figure  3.  Layout  of  the  main  components  of  a  pulse-type  laser  rangefinder  (From 
Shan  &  Toth,  2009) 

Once  one  has  the  components  that  make  up  a  laser,  all  ranging,  profiling,  and 
scanning  are  based  on  laser  ranging  instrument  that  can  measure  distance  to  a  very  high 
degree  of  accuracy.  There  are  two  main  methods  of  measuring  range  using  laser:  the 
time  of  flight  (TOF)  method  and  the  phase  comparison  method.  The  TOF  method  very 
accurately  measures  how  long  it  takes  a  laser  pulse  to  travel  from  the  emitter  to  the  target 
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and  back  to  a  receiver.  The  phase  comparison  method  measures  the  phase  difference  in 
the  sinusoidal  pattern  emitted  by  a  continuous  wave  (CW)  laser.  Due  to  the  limited 
power  of  CW  lasers  there  are  very  few  phase  comparison  systems  used  in  airborne  or 
spaceborne  topographic  mapping  (an  exception  is  the  ScaLARS  research  laser  scanner) 
and  will  not  be  discussed  further  (Shan  &  Toth,  2009). 

The  TOF  method  for  LiDAR  is  analogous  to  range  equations  for  radar  systems 
the  measure  the  precise  travel  time  for  an  emitted  pulse.  The  relationship  between  the 
range  to  the  target  and  the  time  for  the  pulse  to  make  the  round  way  trip  to  the  target  and 
back  to  the  receiver  is: 


R= 


CT 

~1 


Where  R  is  the  slant  range  to  the  target,  c  is  the  speed  of  light,  which  is  known,  and  x  is 
the  measured  time  interval.  Since  the  speed  of  light  is  very  accurately  known,  the 
accuracy  of  the  range  measurement  depends  largely  on  the  precision  of  the  time 
measurement  and  the  method  of  defining  the  leading  edge  of  the  returning  pulse  (no 
rectangular  pulses)  (Baltsavias,  1999).  For  example,  to  achieve  a  1  cm  resolution  the 
timer  should  be  able  to  measure  a  66ps  interval  which  requires  a  clock  rate  of 
approximately  15GHz  (Shan  &  Toth,  2009). 

As  with  RADAR,  technology  approaches  have  devised  to  make  longer  pulses 
more  practical.  The  range  equation  is  a  tool  that  is  widely  used  in  the  literature  that 
relates  the  power  received  (Pr)  after  being  reflected  from  a  target  from  a  laser  pulse  with  a 
given  power  transmitted  (Pt).  Baltsavias  (1999)  provides  an  explanation  of  the  following 
equation: 


P  =  p  ^  A'  P 

'  P  ?rR2  t 


Where  p  is  the  reflectivity  of  the  target,  q  is  atmospheric  transmission,  Ar  is  the  area  of 
the  receiving  optics,  and  R  is  the  range  form  the  LiDAR  sensor  to  the  target.  This 
equation  assumes  the  laser  footprint  and  detector  (optics)  completely  covers  the  target 
area  and  the  reflected  power  radiates  unifonnly  in  a  hemisphere.  This  simplified 
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equation  is  an  optimistic  approximation  since  several  other  factors  internal  to  specific 
LiDAR  systems  introduce  additional  losses,  but  the  equation  does  show  that  at  increased 
sensor  heights  (greater  range)  PTand  Armust  be  increased  to  get  useful  signal. 

As  an  example,  the  Optech  3100  system  that  collected  the  data  for  this  study  has  a 
peak  transmission  power,  PT=  980  watts  and  a  standard  operating  range  of  3,000  ft,  and 
receiver  area  of  10cm  (Optech,  2011).  Assuming  a  typical  mid-latitude  transmittance  of 
0.6  and  reflectivity  of  0.5,  both  for  wavelength  1064  mn  the  above  equation  yields: 


Pr 

Pr 


Q5(0.6)2(0.01m2) 
7t(9  14m)2 
=  6.72x  10~7  W 


980W 


Depending  on  the  distances  being  measured,  it  is  possible  to  make  a  large  number  of 
measurements  over  a  short  period  of  time.  This  is  known  as  the  pulse  repetition 
frequency  and  for  most  commercial  systems  typically  ranges  from  33  to  167  kHz.  To 
continue  with  the  Optech  3100  example:  if  a  distance  of  3,000  ft  (914  m)  needs  to  be 
measured,  then  the  travel  time  over  the  6,000  ft  out  and  back  of  the  pulse  is: 


T= 


2  *R 


T=  2*(914m)  =6.1xlO-S 


2.998x10* 


m 


PRF=- 


1 


6.1x10  5 


=  164kHz 


Thus  the  maximum  PRF  of  the  system  operating  at  3,000  ft  is  164  kHz,  and  generally  the 
higher  the  operating  altitude  the  longer  time  interval  required  between  successive  pulses 
(Shan  &  Toth,  2009).  However,  multiple  pulse  techniques  are  starting  to  be  introduced 
which  allows  for  more  than  one  pulse  to  be  in  the  air  at  a  time,  known  as  multiple  pulses 
in  air  (MPIA).  These  systems’  PRF  are  no  longer  limited  by  altitude  and  pulse  time-of- 
flight,  but  only  limited  by  how  frequently  a  pulse  can  be  emitted  from  the  laser  source. 
As  a  practical  example,  the  Optech  Pegasus  advertises  a  PRF  of  500  kHz  at  an  operating 
altitude  of  300  m  or  225  kHz  at  1,000  m  altitude  (Optech,  2011). 
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3.  Applications  of  LiDAR 

LiDAR  systems  are  generally  used  in  one  of  three  distinct  realms  based  on  the 
platform  of  the  laser  scanner:  terrestrial,  airborne,  and  spaceborne.  Each  type  of 
application  has  distinct  processes,  and  based  on  system  advantages  and  disadvantages  are 
used  by  different  scientific  and  surveying  communities. 

Since  the  beginning  of  laser  range  finders,  surveyors  have  used  light  ranging 
systems  for  distance  and  range  measurements  (Shan  &  Toth,  2009).  In  the  community  of 
field  surveying,  lasers  began  to  replace  tungsten  and  mercury  vapor  lamps  in  electronic 
distance  measuring  instrument  in  the  1970s.  These  instruments  were  only  used  initially 
for  measuring  distances  for  control  surveys  or  geodetic  networks,  but  theodolites  were 
used  separately  to  measure  the  angles  needed  for  these  operations.  Eventually,  these 
separate  systems  were  consolidated,  and  with  the  advent  of  small  eye-safe  lasers, 
reflectorless  distance  measurements  became  possible  (Gordon  &  Charles,  2008).  This 
evolution  in  the  different  field  surveying  applications,  led  to  the  current  tripod  and 
vehicle-mounted  laser  scanning  systems  being  used  for  topographic  mapping 
applications,  such  as  Google  Maps  Street  View,  which  uses  three  lasers  to  capture  3D 
data  (Google). 

For  spaceborne  LiDAR  systems,  the  challenging  operational  environment  has 
made  development  slow.  Spaceborne  LiDAR  systems  function  under  the  same  principles 
described  earlier,  however,  with  distances  and  speeds  100  times  greater  than  that  of  an 
airborne  LiDAR  system.  Much  more  powerful  lasers  are  required  and  system  PRFs  are 
much  reduced.  For  these  reasons,  spaceborne  LiDAR  systems  to  date  have  been  laser 
profiling  systems  and  not  laser  scanning  systems,  only  measuring  pulses  along  track  of 
the  spacecraft  with  no  side-to-side  swath  measurements.  Due  to  pressures  from  various 
scientific  communities  that  were  concerned  with  accurate  elevation  data  for  ice-covered 
terrain  and  desert  regions,  NASA  began  several  Space  Shuttle  experiments  beginning 
with  the  mission,  LITE  (LiDAR  In-space  Technology  Experiment)  (NASA,  1994). 
Conducted  in  1994,  onboard  STS-64,  the  LITE  payload  was  a  2  metric  ton  system  that 
operated  on  three  different  wavelengths:  1064  nm,  532  mn,  and  355  mn  and  was 

principally  concerned  with  atmospheric,  climatic,  and  weather  research.  The  success  of 
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LITE  led  to  two  additional  shuttle  experiments  in  1996  and  1997,  named  Shuttle  Laser 
Altimeter  (SLA).  The  first  mission,  SLA-01,  was  on  STS-72  and  operated  for  over  80 
hours  measuring  over  475,000  pulses  returned  from  land  surfaces.  The  tracks  for  SLA- 
01  are  in  Figure  4. 
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Figure  4.  Ground  track  of  STS-72  while  collecting  data  for  SLA-0 1 .  Color  ramp 
indicates  altitude  in  meters  (From  Harding,  2001) 

The  second  experiment  took  place  on  STS-85  in  August,  1997.  SLA-02 
generated  100  m  laser  footprints  on  the  ground  at  a  PRF  of  10  Hz  and  wavelength  of 
1064nm  giving  a  ground  spot  spacing  of  700  m  (Carabajal  et  ah,  1999).  These 
experiments  were  the  precursors  to  the  Geoscience  Laser  Altimeter  System  (GLAS)  that 
is  the  primary  payload  on  the  Ice,  Cloud,  and  land  Elevation  Satellite  (ICESat). 
Launched  in  January  2003,  GLAS  has  the  primary  mission  of  monitoring  the  Earth’s  ice 
sheets  to  determine  changes  in  total  mass  and  any  contribution  to  sea  level  changes 
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(Schutz,  Zwally,  Shuman,  Hancock,  &  DiMarzio,  2005).  GLAS  operates  at  two 
wavelengths,  1064  nm  and  532  nm,  and  from  an  altitude  of  600  km  creates  a  70  m  laser 
spot  along  the  ground  track  at  40  Hz  PRF  resulting  in  170  m  interval  pulses  (Abshire  et 
ah,  2005).  ICESat  was  decommissioned  on  August  17,  2010,  after  a  successful  seven 
years  of  operations.  Figure  5  describes  tree  canopy  height  in  the  United  States  from  data 
collected  by  GLAS. 


Figure  5.  A  forest  canopy  height  map  of  the  contiguous  United  States  (From  Lefsky, 
2010) 

Airborne  laser  scanning  in  recent  years  has  developed  into  a  powerful  tool  used 

with  much  success  in  the  fields  of  engineering,  archeology,  and  forestry  (Vosselman  & 

Maas,  2010).  In  these  fields,  airborne  LiDAR  systems  have  shown  some  clear 

advantages  to  traditional  data  collection  methods.  LiDAR  has  the  capacity  for  high  data 
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densities  and  good  data  accuracy.  Some  of  the  highest  measurement  densities  are 
captured  by  helicopter  mounted  systems  due  to  their  slower  airspeed  and  can  take  up  to 
30  measurements/m  .  Fast  data  acquisition  is  another  advantage  of  airborne  systems,  as 
opposed  to  field  measurements  on  the  ground.  Airborne  LiDAR  systems  also  require  a 
minimal  amount  of  ground  truth.  Field  work  is  minimal  since  LiDAR  generally  only 
requires  a  few  ground  reference  points,  even  for  large  data  collection  areas  (Vossehnan  & 
Maas,  2010). 

For  engineering  practices,  airborne  LiDAR  systems  are  used  in  corridor  mapping 
of  outdoor  structures  such  as  railroads,  power  lines,  pipelines,  and  dikes.  For  power  line 
and  railroad  monitoring,  LiDAR  is  useful  to  map  and  then  model  the  immediate  area 
surrounding  the  subject  of  interest.  Vegetation  growing  too  close  can  be  identified  and 
removed.  For  pipeline  studies,  terrain  around  the  pipelines  is  studied  to  assess  the  impact 
and  possible  damage  that  could  occur  to  the  pipeline  from  the  terrain.  Dike  and  levee 
monitoring  programs  also  use  LiDAR  to  monitor  and  study  the  structure  profiles  as  well 
as  model  the  effects  from  potential  failures.  For  example,  in  the  aftennath  of  the  collapse 
of  two  minor  dikes  in  the  Netherlands,  the  local  Water  Boards  have  used  helicopter 
LiDAR  systems,  FLI-MAP,  to  regularly  monitor  over  6,500  kilometers  of  embankment 
(Franken  &  Flos,  2005). 

In  Detection  and  Vectorization  of  Roads  from  LiDAR  Data,  (Clode  et  ah,  2007), 
road  classification  is  conducted  for  urban  terrain  to  support  city  planning  and  road 
network  updating.  They  propose  a  method  for  automatic  detection  and  vectorization  of 
road  networks  solely  from  LiDAR  data.  For  the  purposes  of  their  research,  roads  are 
assumed  to  be  on  the  digital  terrain  model  (DTM);  that  is,  only  the  LiDAR  points  within 
a  given  tolerance  of  the  DTM  are  considered  as  candidates  for  road.  Next,  a  training 
algorithm,  based  on  road  intensity  values,  is  used  further  classify  roads  in  the  scene. 
Morphological  filtering  is  then  conducted  to  remove  gaps  caused  by  overhanging  trees  or 
reflections  from  objects  such  as  vehicles  on  the  roads.  At  this  point,  the  classified  road 
images  are  vectorized  using  a  convolution  of  the  binary  image  with  a  phase  coded  disk 
(PCD).  The  road  centerline,  orientation,  and  width  are  determined  by  convolution  with 
the  PCD  and  then  vectorized.  Typical  results  are  shown  in  Figure  6. 
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Figure  6.  Fairfield,  New  South  Wales,  road  centerline  (left)  and  edges  (right)  are 

overlaid  on  aerial  imagery  to  demonstrate  results  (From  Clode  et  ah,  2007) 

Airborne  LiDAR  shows  great  promise  as  a  tool  to  inspect  archeological  sites  in 
forested  areas.  Woodlands  can  protect  archeological  remains  from  erosion,  but  also 
conceal  archeological  structures  that  could  otherwise  be  identified  from  aerial 
photography.  In  these  wooded  areas,  LiDAR  can  create  detailed  DTMs  that  have  made 
important  contributions  to  archeological  prospection  (Vossehnan  &  Maas,  2010).  The 
uses  of  LiDAR  in  archeological  survey  have  some  of  the  same  considerations  as  using 
LiDAR  to  identify  trails  under  tree  canopy.  To  identify  features  beneath  tree  canopy,  a 
high  initial  point  density  is  required  to  ensure  a  sufficient  number  of  pulse  measurements 
are  returned  from  the  ground.  In  an  archeological  study  of  the  Mayan  City  of  Caracol  the 
LiDAR  collection,  depicted  in  the  figure  below,  consisted  of  two  sets  of  perpendicular 
flight  lines  with  50%  overlap,  low  operating  altitude  (800  m),  and  a  high  PRF  (100  kHz) 
which  resulted  in  an  average  of  20  returns  per  square  meter  (Diaz,  2011).  This  area  has 
been  heavily  ground  truthed  and  studied  using  traditional  archeological  survey  techniques 
by  researchers  at  the  University  of  Central  Florida  since  1983.  After  comparisons  of  the 
LiDAR  data  to  the  ground  surveys,  the  agreement  was  very  good,  and  over  the  study  area 
the  LiDAR  data  showed  higher  spatial  resolution  identifying  many  new  archeological 
features  previously  undiscovered  by  ground  archeology  teams  (Diaz,  2011). 


Figure  7.  Map  showing  the  Vaca  Plateau  project  area  in  western  Belize  defined  by  the 
white  box.  The  colored  lines  represent  the  ground  tracks  of  the  different 
flights  (From  Diaz,  2011) 

An  elegant  treatment  of  LiDAR  use  in  archeology  was  written  by  Crutchley  and 
Crow  (2010)  for  English  Heritage,  the  United  Kingdom’s  statutory  advisor  on  the  historic 
environment.  They  discuss  LiDAR  fundamentals,  file  format  types,  project  planning, 
data  manipulation,  and  interpretation.  Of  particular  interest  are  the  five  case  studies 
conducted  and  discussed.  The  case  studies,  which  were  conducted  from  2005  through 
2008,  built  upon  each  other  going  from  a  previously  intensely  studied  area,  Stonehenge, 
to  a  later  project,  Savernake  Forest,  which  assessed  the  value  of  using  LiDAR  in  a 
woodland  environment  (Crutchley  &  Crow,  2009).  They  found  that  by  using  bare  earth 
hill  shaded  products  with  the  maximum  vegetation  removed;  previously  unrecorded 
archeological  monuments  were  identified  for  further  ground  study. 
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Figure  8.  Enclosures  at  Church  Walk  from  LiDAR  (left)  and  ground  survey  (right) 
revealing  almost  exact  agreement  on  size,  location,  and  shape  of 
archeological  features  (From  Crutchley  &  Crow,  2009) 


Since  airborne  LiDAR  became  commercially  viable  in  the  1990s,  foresters  have 
used  scanning  systems  and  processes  in  the  areas  of  forest  management,  inventory, 
carbon  sink  analysis,  biodiversity  characterization,  and  habitat  analysis.  Conventional 
means  of  conducting  these  studies  involve  physically  walking  the  ground  and  taking  key 
measurements;  tree  height,  diameter,  and  volume,  for  a  small  representative  plot  and  then 
extrapolating  those  results  to  include  the  entire  stand  of  trees.  These  techniques  have  an 
estimated  accuracy  of  only  15-30%  (Vosselman  &  Maas,  2010).  However,  using  LiDAR 
these  same  measurements  can  be  accomplished  to  measure  the  same  parameters  for  full- 
field  rather  than  plot  based  forest  inventories.  Also,  with  increased  point  density 
available  in  modern  systems,  LiDAR  can  conduct  measurements  over  a  whole  tree  stand 
(requiring  point  spacing  of  2-4  m)  or  can  make  measurements  on  a  tree-wise  approach 
for  individual  trees  (requires  more  than  one  point  per  square  meter). 

Within  the  field  of  forestry,  forest  road  planning  and  mapping  plays  an  important 
role  in  management  activities.  Forest  management  professionals  can  use  road  inventory 
data  to  evaluate  land  use  impacts,  watershed  disturbances,  and  future  planning  of  forest 
road  networks  (White,  2010).  Electro-optical  satellite  imagery  and  likewise  aerial 
photography  have  not  been  effective  methods  for  forest  road  mapping  since  even  sparse 
tree  cover  prevents  those  sensors  from  observing  the  trails  under  canopy.  Ground 
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inventories  of  forest  roads  are  time  consuming,  expensive,  and  subject  to  much 
uncertainty.  Across  studies  and  estimates  of  forest  service  road  systems,  from  20%  up  to 
as  much  as  50%  of  road  lengths  are  thought  missing  from  inventories,  and  roads  and 
trails  in  privately  held  lands  are  seldom  inventoried  or  reported  (R.  A.  White,  Dietterick, 
Mastin,  &  Strohman,  2010).  In  their  research,  White  et  al.  (2010)  manually  mapped  a 
forest  road  beneath  tree  canopy  to  a  positional  accuracy  of  1.5  m  and  road  grade 
measurements  within  0.53%  of  ground  survey.  Similarly,  Espinoza  and  Owens  (2007) 
demonstrated  an  82%  accuracy  rate  of  covered  trail  identification  from  LiDAR  derived 
products. 


Figure  9.  Field  surveyed  centerline  (black)  and  digitized  centerline  (red)  (From  White 
et  al.,  2010) 

David  et  al.,  2009,  attempt  to  automate  the  forest  road  detection  process  using 
only  FiDAR  data.  Following  a  similar  methodology  as  (Clode  et  al.,  2007),  the 
researchers  use  three  FiDAR  derived  products  to  detect  trails.  The  three  products  are  a 
nonnalized  digital  surface  model  (nDSM),  an  altimetric  variance  image,  and  an  intensity 
image.  Using  these  three  layers  the  researchers  identified  mean  and  variance  for  trails  in 
each  or  the  layers,  and  then  used  a  region  growing  algorithm  to  expand  trail  sections  from 
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a  researcher  selected  seed.  The  resulting  binary  trail  image  is  vectorized  using 
morphological  tools,  and  finally  the  trail  borders  and  centerlines  are  exported  as  vector 
data.  Results  from  their  study  are  shown  in  Figure  10. 


Figure  10.  Results  of  pathway  vectorization.  Red  is  detected  and  vectorized  pathways, 
pink  are  detected  borderline  points,  and  blue  are  pathways  from  an  existing 
ground  surveyed  database  (From  David  et  ah,  2009) 


B.  THEORY 

1.  Supervised  Classification  Techniques 

The  purpose  of  this  research  is  to  identify  roads  and  trails  under  canopy.  The 
advantage  of  LiDAR  is  its  ability  to  “see  to  the  ground”  and  derive  digital  terrain  model 
(DTM)  products  to  analyze.  Once  these  products  are  created  the  problem  of  trail 
detection  becomes  an  image  processing  problem,  and  several  techniques,  traditionally 
used  on  multispectral  and  hyperspectral  images,  can  be  applied  to  the  LiDAR  derived 
products.  Supervised  classification  generally  refers  to  a  set  of  algorithms  used  to  classify 
an  image  based  on  statistical  analysis.  Regardless  of  the  statistical  methods  used,  all 
classification  techniques  generally  follow  the  same  steps.  First,  decide  on  the  number  of 
classes  the  image  is  to  be  classified  into,  for  this  research  the  number  of  classes  is  two: 
trail  and  not  trail.  Second,  pick  a  set  of  pixels  in  the  image  that  represent  the  typical 
characteristics  of  the  classes  from  step  one.  These  pixels  are  known  as  the  training  data 
for  the  classifier  (Richards  &  Jia,  2006).  Third,  this  training  data  is  used  to  determine  the 
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parameters  used  for  a  probability  model  or  some  other  set  of  equations  that  quantitatively 
describe  the  classes  from  step  one.  Next,  using  the  trained  classifier,  label  every  pixel  in 
the  image  as  belonging  to  one  of  the  classes  from  step  one,  from  this  produce  class  maps 
showing  the  results  of  the  classification  or  produce  the  results  in  tabular  fonn.  Finally, 
compare  the  results  of  the  class  maps  to  ground  truth  obtained  through  site  visits,  survey 
data,  of  even  by  photo  interpretation  of  the  images.  This  process  can  be  reiterated  to 
improve  classification  results  by  refining  the  training  data  parameters  based  on 
knowledge  of  the  original  classes. 

2.  Maximum  Likelihood  Classification 

Maximum  likelihood  classification  is  the  most  common  type  of  supervised 
classification  used  with  remote  sensing  images  (Richards  &  Jia,  2006).  This 
classification  method  is  based  on  mean  statistics;  this  method  uses  a  Bayesian  probability 
function  that  has  been  calculated  from  the  training  data  classes.  Then  each  pixel  is 
judged  as  to  which  class  it  most  probably  belongs.  The  following  approach  is  developed 
and  explained  in  sufficient  detail  for  most  remote  sensing  applications  by  Richards 
(2006),  but  the  process  comes  from  the  field  of  mathematical  pattern  recognition  and 
machine  learning,  and  is  covered  in  greater  mathematical  detail  in  that  discipline. 

For  each  of  the  classes  in  the  training  data  in  the  image 

&>.,/  =  1  ,...M 

where  M  is  the  number  of  classes.  To  detennine  to  which  class  a  pixel  vector  x  belongs 
the  conditional  probabilities  are  used. 

p{coi  | jc), /  =  1  ,..M 

Vector  x  is  a  measurement  of  the  brightness  of  a  given  pixel  describing  its  location  in 
multispectral  space.  The  probability  p(coi  |jc)  is  the  likelihood  that  the  class  (>)I  i s  correct 
for  any  given  pixel  at  jc  .  At  this  point  classification  can  be  performed  by 

x  e  if  p(ooi  |jc)  >  p(a>j  |jc)  for  all  j  ^  i 
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stating  that  a  pixel  at  x  belongs  to  class  op  if  p{coi  | _v)  is  the  largest.  However,  p(mi  jc) 
is  unknown  and  must  be  computed  by  Bayes’  theorem: 

p(op  |x)  =  p{x\a>i)p{(ot)l  p(x) 

where  p(x |o.)  is  estimated  from  the  pixels  in  the  training  data,  p(op)  is  the  probability 
that  <z>.  is  in  the  image,  and  p(x)  is  the  probability  of  finding  a  pixel  from  any  class  at 
location  x ,  and  can  be  computed  by 

M 

p(x)  =  YJP(x\a,)p(o)i) 

i= 1 

Since  p{x\a>j)  are  known  from  the  training  data  and  p{(oi )  may  be  estimated  from 
knowledge  of  the  image  it  is  more  acceptable  to  perfonn  the  classification  as  follows: 

jceffl;,  if  p{x\a>j) p{a>i)>  p{x\a>j) p{a>j)  for  all  j^i 

In  this  expression  p(x)  has  been  removed  as  a  common  expression. 

To  summarize,  assuming  the  pixels  in  each  class  are  normally  distributed  in 
multidimensional  space,  the  maximum  likelihood  classifier  computed  both  the  variances 
and  covariance  of  the  training  data  when  assigning  pixels  to  one  of  the  classes  in  the 
training  data.  Most  image  processing  software  allows  the  analyst  to  assign  training  sets 
and  performs  the  above  classification  in  the  multidimensional  space  of  the  image.  The 
above  treatment  will  assign  every  pixel  to  one  of  the  specified  classes,  however,  if  not  all 
unique  classes  in  an  image  were  identified  in  the  training  data  this  may  lead  to  a  poor 
classification  result  (Richards  &  Jia,  2006).  For  this  reason,  many  software  packages  also 
allow  the  analyst  to  set  a  threshold  where  any  pixel  below  the  threshold  probability  for  all 
classes  are  not  classified  as  depicted  in  Figure  11. 
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Figure  11.  a  Illustration  of  poor  classification  for  patterns  lying  near  the  tails  of  the 
distribution  functions  for  all  of  the  classes;  b  Use  of  thresholds  to  remove 
poorly  classified  regions  (From  Richards  &  Jia,  2006) 

3.  Edge  Detection  and  Enhancement 

Another  method  explored  in  the  course  of  this  research  was  the  use  of 
neighborhood  operations  to  enhance  the  trails  in  the  LiDAR  derived  products.  These 
procedures  modify  the  brightness  of  an  image  pixel  as  a  function  of  some  weighted 
average  of  the  brightness  of  the  surrounding  pixels  (Richards  &  Jia,  2006).  The  methods 
attempted  were  template  operators,  where  a  box  or  window  size  is  defined  (template)  and 
then  moved  over  the  original  image  row  by  row  and  column  by  column.  The  template  is 
also  known  as  the  kernel.  The  operation  takes  the  product  of  the  pixel  values  of  the 
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portion  of  the  original  image  covered  by  the  template  and  the  template  values,  and  these 
values  are  summed.  The  summed  value  is  the  new  brightness  value  for  the  center  pixel 
of  the  original  image. 


45 

60 

98 

127 

132 

133 

137 

133 

46 

65 

98 

123 

126 

128 

131 

133 

47 

65 

96 

115 

119 

123 

135 

137 

47 

63 

91 

107 

113 

122 

138 

134 

50 

59 

so 

97 

110 

123 

133 

134 

49 

53 

68 

83 

97 

113 

128 

133 

50 

50 

58 

70 

84 

102 

116 

126 

50 

50 

52 

58 

69 

86 

101 

120 

69 

95 

116 

125 

129 

132 

68 

92 

110 

120 

126 

132 

66 

86 

104 

114 

124 

132 

62 

78 

94 

108 

120 

129 

57 

69 

83 

98 

112 

124 

53 

60 

71 

85 

100 

114 

0.1 

0  1 

0  1 

0.1 

0.2 

0.1 

0.1 

0  1 

0.1 

Figure  12.  Neighborhood  filtering  (convolution):  the  image  on  the  left  is  convolved 
with  the  filter  in  the  middle  to  yield  the  image  on  the  right.  The  light  blue 
pixels  indicate  the  source  neighborhood  for  the  light  green  destination  pixel 
(After  Szeliski,  2011) 


These  templates  can  be  used  to  filter  an  image,  adding  blur,  sharpening, 
enhancing  edges,  or  removing  noise  in  either  color  images  or  binary  images  as  seen  in 
Figure  13.  For  this  research,  several  filters  were  experimented  with,  however,  a 
Laplacian  convolution  filter  was  found  to  be  most  effective.  It  is  a  second  derivative 
edge  enhancement  filter  that  operates  without  regard  to  the  direction  of  the  linear  feature 
being  enhanced.  It  uses  a  kernel  with  a  high  central  value,  typically  negative  values  in 
the  up-down  and  side-to-side  directions,  and  values  of  zero  on  the  corners  of  the  filter. 
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Figure  13.  Examples  of  neighborhood  operations:  (a)  original  image;  (b)  blurred;  (c) 
sharpened:  (d)  smoothed  with  edge-preserving  fdter;  (e)  binary  image;  (f) 
dilated;  (g)  distance  transform;  (h)  connected  components.  For  dilation  and 
connected  components,  black  pixels  are  assumed  to  have  a  value  of  1  (From 
Szeliski,  2011) 


4.  Decision  Trees 

The  above  techniques  are  considered  single  stage  classifications  in  that  only  one 
decision  or  operation  is  performed  on  a  pixel.  An  alternate  method  is  using  a  decision 
tree,  where  a  series  of  binary  decisions  are  made  to  determine  the  correct  category  for 
each  pixel  (Richards  &  Jia,  2006).  The  advantage  to  using  the  decision  tree  approach  is 
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that  at  each  level  of  the  decision  tree  different  data  sets,  different  class  attributes,  and 
even  different  algorithms  can  be  used  to  improve  the  accuracy  of  the  overall 
classification  of  each  pixel. 

Using  the  decision  tree  process  both  the  maximum  likelihood  classification  and 
the  edge  enhancement  can  be  performed  on  the  same  data  set.  Using  a  large  to  small 
approach,  the  maximum  likelihood  classification  identifies  all  pixels  that  have  the  same 
characteristics  are  the  trail  training  set.  These  trail  pixels  are  labeled  as  the  survivors  of 
node  one  of  the  tree,  those  survivors  move  on  to  the  next  node  of  the  decision  tree  where 
the  previously  described  Laplacian  edge  enhancement  filter  can  be  applied  only  to  the 
survivors  mask  from  node  one.  The  resulting  survivors  from  node  two  have  been 
identified  as  trail  by  both  the  maximum  likelihood  classification  and  the  edge 
enhancement.  At  this  point  any  erosion  or  dilation  operators  can  be  applied  to  the 
survivors  from  node  two  in  order  to  reduce  noise  or  enhance  linear  features,  respectively. 
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III.  PROBLEM 


A.  OVERVIEW 

The  location  of  the  LiDAR  data  set  for  this  study  is  in  the  Santa  Cruz  Mountains 
in  California.  Swanton  Pacific  Ranch  is  19  km  north  of  Santa  Cruz  and  is  owned  and 
managed  by  California  Polytechnic  State  University  Corporation  under  the  College  of 
Agriculture,  Food,  and  Environmental  Sciences  as  an  educational  and  research  facility 
Swanton  Pacific  Ranch  is  a  subset  of  the  Scotts  Creek  Watershed,  and  the  southern 
portion  of  the  Little  Creek  Watershed  is  contained  in  Swanton  Pacific  Ranch.  The 
topography  of  the  study  area  is  varied  from  the  lower  part  of  the  property  that  is  near  sea 
level  to  the  hills  of  the  ranch  that  reach  elevations  of  488  m.  The  extent  of  the  study  area 
has  an  average  slope  of  45%. 
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Figure  14.  The  Little  Creek  watershed,  tributary  to  the  Scotts  Creek  watershed. 

Swanton  Pacific  Ranch  property  boundary  in  red  (From  White,  2010) 


The  forest  canopy  in  the  study  area  consists  primarily  of  second-growth  coast 
redwood,  but  also  has  Douglas-fir  and  mature  red  alder.  The  overstory  of  the  study  has 
been  measured  using  a  vertical  densitometer  at  thirty  forest  inventory  locations 
throughout  Swanton  Pacific  Ranch  to  support  forest  management.  The  canopy  coverage 
at  these  sites  range  between  40  to  96%  density,  averaging  80%  (White,  2010). 
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Figure  15.  Overstory  examples  in  Swanton  Pacific  Ranch  over  a  forest  haul  road  on 
left  and  over  a  foot  trail  on  right 


Typical  roads  and  trails  in  the  area  range  from  forest  hauls  roads  that  follow 
fonner  railroad  grades  to  smaller  foot  trails  that  are  approximately  one  meter  wide. 
Examples  of  typical  roads  and  trails  identified  on  hillshade  products  are  shown  in  Figure 
16. 
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Figure  16.  Typical  road  on  left  with  2  meter  tape  measure  and  typical  trail  on  right 
with  1  meter  tape  measure 

B.  DATA  SET  AND  COLLECTION  METHOD 

The  primary  data  set  for  this  research  was  contracted  by  Cal  Poly,  and  flown  by 
the  Airborne  1  Corporation  in  March  of  2010.  The  LiDAR  system  used  was  an  Optech 
3100/(EA)  operating  at  3,000  ft  altitude  above  ground  level  (AGL)  with  a  PRF  of  150 
kHz.  The  LiDAR  survey  parameters  are  summarized  in  Table  1. 
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Table  1 .  S wanton  Pacific  Ranch  LiDAR  collection  parameters 

(From  Airborne  1,  2010) 


LiDAR  Survey  Parameters 

Aircraft: 

Navajo  Chieftain 

Sensor: 

Op  tech  3100/(EA) 

Altitude: 

3,000  ft  AGL 

Scan  angle: 

14  degrees 

Scan  Frequency: 

30  Hz 

Pulse  repetition  frequency: 

150  kHz 

Returns: 

4  per  pulse  (with  intensity) 

The  extent  of  the  LiDAR  data  covers  approximately  11,200  acres,  encompassing 
all  of  the  Little  Creek  watershed,  and  Swanton  Pacific  Ranch  property  east  of  Scotts 
Creek,  as  shown  in  Figure  15.  The  data  was  delivered  in  the  American  Society  for 
Photogrammetry  and  Remote  Sensing  (ASPRS)  standardized  format  of  LAS  1.1  files. 
The  data  set  consists  of  248  separate  LAS  files  ranging  in  file  size  from  852  KB  to  247 
MB. 

For  this  project  the  vendor,  Airborne  1  Corporation,  supplied  all  LAS  files  in  the 
California  State  Plane  Zone  3,  NAD83  for  the  XY  extent  and  NAVD88  for  the  Z  extent, 
in  U.S.  survey  feet.  The  advantage  of  using  this  is  due  to  the  relatively  small  area  of 
study,  the  distortions  in  the  state  plane  coordinate  systems  are  less  than  Universal 
Transverse  Mercator  (UTM)  or  other  common  projections.  However,  through  the  course 
of  the  project  and  while  working  with  multiple  software  suites  to  analyze  the  same 
dataset,  some  smaller  analysis  packages  do  not  support  state  plane  projections.  To  work 
around  this,  the  data  is  analyzed  in  the  smaller  software  package,  for  example  QT 
Modeler  by  Applied  Imagery,  and  those  analysis  results  are  then  exported  as  either 
geotiffs  or  rasters,  which  can  then  be  reprojected. 
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Figure  17.  Swanton  Pacific  Ranch  boundary  in  red  and  extent  of  LiDAR  data  coverage 
in  blue  (Imagery  from  Google  Earth) 

C.  SOFTWARE  USED 

1.  Economic  and  Social  Research  Institute  (ESRI)  ArcGIS 

ArcGIS  Desktop  is  a  software  product  created  by  ESRI.  It  is  used  by  many  GIS 
professionals  to  compile,  use,  and  manage  geographic  infonnation  in  both  a  standalone 
version  and  for  server  enterprises  (ESRI).  ArcMap  is  one  of  the  desktop  products  within 
the  ArcGIS  suite,  and  is  used  to  construct  maps,  spatial  analysis,  and  data  compilation. 
The  Swanton  Pacific  Ranch  LiDAR  data  set  was  used  with  ArcMap  version  10.0  to 
create  a  consolidated  DEM  from  the  248  LAS  files. 

2.  Environment  for  Visualizing  Images  (ENVI) 

There  are  a  number  of  different  software  suites  for  processing  geospatial  imagery. 
The  classification  and  morphological  operations  for  this  thesis  were  conducted  using 
ENVI  from  ITT  Visual  Information  Solutions.  ENVIs  built-in  topographic  modeling, 
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classification,  decision  tree  builder,  and  morphological  operations  were  used  to  conduct 
analysis.  Further  explanation  of  the  use  of  these  tools  is  provided  in  the  Methods  and 
Observations  Chapter. 

D.  FIELD  EQUIPMENT 

Table  2  summarizes  the  additional  equipment  used  to  collect  ground  truth  data. 


Table  2.  Field  equipment 


Field  Equipment 

Equipment 

Description 

Gannin  GPSMAP  60CSX 

Hand-held  GPS  receiver  to  collect  ground  track 
infonnation  and  verify  control  point 

Trimble  Nomad 

Outdoor  rugged  hand-held  computer/GPS  receiver  to 
collect  control  points 

Olympus  Stylus  1030  SW 

Digital  camera  to  capture  trail  characteristics  and 
overhead  cover 
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IV.  METHODS  AND  OBSERVATIONS 


A.  DEM  CREATION 

The  data  set  and  point  density  were  carefully  evaluated  for  this  research  to  ensure 
accurate  results  from  products  derived  from  the  LiDAR  point  cloud,  specifically  ground 
point  density  as  it  affects  the  accuracy  for  DEMs  and  DSMs  particularly  under  forest 
canopy.  In  ArcMap,  the  Point  File  Information  tool  under  the  3D  Analyst  toolbox 
computes  point  spacing  for  each  LAS  file  by  comparing  the  size  of  the  XY  extent  of  the 
file  to  the  file’s  point  count.  This  tool  works  to  quickly  summarize  large  data  sets  since  it 
can  obtain  all  needed  information  from  the  LAS  header.  As  seen  in  Figure  18,  the  data 
set  for  all  returns  ranges  from  0.55  to  1.77  ft  between  points,  with  a  mean  of  0.77  ft. 
However,  if  only  ground  returns  are  considered,  the  point  spacing  increases  to  a  range  of 
1.13  to  4.62  ft,  with  a  mean  of  1.99  ft. 


Swanton  Pacific  Ranch 
LiDAR  Point  Spacing 
All  Returns 


NAD83.  C  lAfomiv  SWI*  Plane  Zone  If 
Denved  from  enaaery  provided  by  Remole 
Senyng  Cenlii.Noat  Postgraduate  SdiOni 


Legend 


Point  Spacing  (ft) 

1 1  13- 1  54 
1 1  55  •  1  94 
U  1  95  ■  2 44 
|?45  -  3  ?1 
1 3  22  •  4  62 


Swanton  Pacific  Ranch 
LiDAR  Point  Spacing 
Ground  Returns 


UAD83.  California  Stale  Plane  Zune  III 
Dem»d  from  ro*aerjr  provided  Remote 
Sensing  Center,  Nmal  Poslgraduele  Sctiool 


Figure  18.  Point  density  comparison  between  all  returns  and  only  ground  classified 
returns 
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At  this  point,  sufficient  point  density  exists  to  begin  DTM  creation.  For  DTM 
creation  in  ArcMap,  the  workflow  begins  by  organizing  the  data  in  the  3D  point  cloud 
into  a  format  that  ArcMap  can  read  (ESRI,  2010).  To  do  this,  convert  the  LAS  files  to  a 
Multipoint  Feature  set  using  the  LAS  to  Multipoint  tool  found  in  3D  Analyst  -> 
Conversion  ->  From  File.  To  build  the  terrain  dataset  only  the  ground  classified  points 
are  used  for  the  terrain  wizard.  The  terrain  wizard  is  found  under  a  File  Geodatabase 
Feature  Dataset,  by  right  clicking  and  choosing  New  Terrain.  To  create  the  DEM  use 
the  Terrain  to  Raster  tool  found  in  3D  Analyst  ->  Conversion  ->  From  Terrain;  creates 
the  output  DTM  using  the  following  parameters:  Method:  Natural_Neighbors,  Cellsize: 
3.28,  and  Pyramid  Level:  0. 


Figure  19.  Terrain  to  Raster  tool  in  ArcMap  used  for  DTM  creation  of  1  meter 
resolution 

Method  refers  to  method  of  interpolation,  and  choices  are  natural  neighbor  or 
linear.  These  methods  are  TIN-based  interpolation  methods  applied  through  the 
triangulated  terrain  surface.  Natural  neighbor,  while  not  as  fast  as  linear,  generally 
produces  better  results  in  terms  of  accuracy.  This  method  finds  the  closest  set  of  input 
points  to  a  query  point  and  applies  weights  to  them  based  on  proportionate  areas  to 
interpolate  a  value  (Sibson,  1981).  The  cell  size  sets  the  output  resolution;  for  this 
example  3.28  produces  a  1  meter  resolution  DTM  since  the  original  data  set  is  in  U.S. 
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feet.  The  resolution  parameter  indicates  which  pyramid  level  of  the  terrain  dataset  to  use 
for  conversion.  To  output  a  raster  dataset  at  full  resolution,  this  parameter  is  set  to  0. 


Swanton  Pacific  Ranch 
Digital  Terrain  Model  (DTM) 


MAD83.  California  Slal*  Plan*  Zone  III 
Oem<n5  licirn  tnnqvry  pnwidpti  by  R«no1« 
Sensing  CmiIM,  Nar>al  Postgraduate  School 


Figure  20.  DTM  at  1  meter  resolution  produced  from  LiDAR  point  cloud 

B.  CLASSIFICATION 

The  1  meter  DTM  was  imported  into  ENVI  for  characterization  and  classification 
of  the  data,  specifically,  defining  trail  parameters  that  could  be  observed  from  the  DTM 
and  DTM  derived  products.  Using  the  ENVI  Topographic  Modeling  tool  from  the  main 
menu,  the  software  computes  statistics  from  the  DTM  and  outputs  1 1  bands,  each 
depicting  a  separate  characteristic  derived  from  the  DTM  as  seen  in  Figure  21. 
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ENVI  Topographic  Modeling  Bands 


Profile  Convexity  Plan  Convexity  Longitudinal  Convexity  Cross  Sectional  Convexity 


Minimum  Curvature  Maximum  Curvature  RMS  Error  Slope  Percentage 


Figure  2 1 .  Full  scene  images  representing  the  eleven  bands  created  by  the  ENVI 
Topographic  Modeling  tool  from  a  DTM 
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The  topographic  modeling  tool  in  ENVI  takes  a  DTM  input  and  extracts 
parameters  of  the  scene.  ENVI  calculates  these  parameters  by  fitting  a  quadratic  surface 
to  the  DTM  for  a  user  entered  kernel  size  and  then  taking  the  appropriate  derivatives  (ITT 
Visual  Information  Solutions,  2010).  The  kernel  size  of  three  was  used  for  this  research. 
The  slope  is  measured  both  in  degrees  and  percentage,  with  zero  degrees  corresponding 
to  a  horizontal  surface,  and  slope  percentage  calculated  by  the  formula:  100*rise/run. 
The  aspect  is  the  direction  that  a  surface  faces,  and  is  calculated  with  zero  degrees  to  the 
north  and  increasing  clockwise.  Several  convexities  are  detennined,  and  the  profile  and 
plan  convexity  measure  the  rate  of  change  of  the  slope  and  aspect  respectively.  The 
longitudinal  convexity  is  the  measure  of  the  surface  curvature  in  the  down  slope 
direction,  and  the  cross-sectional  convexity  is  the  measure  of  the  surface  curvature  in  the 
across-slope  direction.  The  minimum  and  maximum  curvatures  are  calculations  for  the 
overall  scene  surface,  and  the  root  mean  square  (RMS)  error  measure  how  well  the 
quadratic  surface  fits  the  DTM  data. 

At  this  time,  training  data  is  identified  in  the  scene  that  represents  a  typical 
section  of  trail  under  tree  canopy.  Within  the  scene  a  region  of  interest  (ROI)  is 
identified  that  includes  pixels  of  trail  and  nontrail  in  the  scene. 
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Figure  22.  Initial  ROIs  (in  white)  to  use  as  inputs  for  training  data.  Image  depicts 
slope  in  degrees 

The  pixels  contained  in  the  ROIs  were  exported  to  the  n-D  visualizer.  The  n-D 
visualizer  is  used  in  spectral  image  analysis  to  locate,  identify,  and  cluster  the  purest 
pixels  and  extreme  spectral  responses  within  a  data  set.  For  this  research  the  n-D 
visualizer  was  used  to  check  the  separability  of  pixel  clusters  by  visualizing  the  data 
cloud  in  n-D  space  using  the  image  bands  as  plot  axes.  The  data  can  be  rotated,  grouped 
into  classes,  and  those  classes  can  be  collapsed  to  make  finer  classification  easier. 
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Figure  23.  Initial  ROI  pixels  plotted  in  n-D  visualizer  on  left  plotted  in  all  eleven  band 
axes,  user  clustered  pixels  in  n-D  visualizer  on  right  plotted  using  bands  1 , 
2,  and  4 


The  results  of  the  used  clustered  pixel  classes  were  then  exported  back  to  the 
image  as  new  class  ROIs.  In  this  example,  the  green  pixels  represent  training  data  pixels 
for  trail;  all  other  classes  are  considered  not  trail  in  red,  blue,  and  yellow  pixels. 


Figure  24.  Clustered  class  ROIs  overlaid  on  slope  image,  green  pixels  are  trail  and  all 
other  classes  are  not  trail 
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Using  201  green  pixels  as  the  exemplar  for  a  trail  under  canopy,  out  of  a  total  of 
over  14  million  pixels  in  the  scene,  supervised  classification  was  conducted.  Using 
maximum  likelihood  classification,  the  mean  and  variance  for  each  class  in  Figure  24  is 
computed  for  each  of  the  eleven  bands.  Then,  each  pixel  is  classified  into  one  of  the  four 
classes  to  which  it  has  the  highest  probability  of  belonging. 


Figure  25.  Result  from  maximum  likelihood  classification  with  all  classes  (upper  left), 
same  spatial  subset  with  only  trails  identified  in  green  (upper  right),  same 
spatial  subset  showing  the  slope  in  degrees  band(lower  left)  and  overhead 
imagery  for  comparison  (lower  right) 
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Rule  images  were  generated  during  the  classification  process  as  well.  The  rule 
images  show  the  classification  results,  before  the  final  assignment  of  classes.  For 
example,  pixels  in  the  rule  image  for  n_D  class  #2  (Trail)  for  the  maximum  likelihood 
classifier  represent  the  likelihood  that  it  belongs  to  that  class.  A  rule  image  for  maximum 
likelihood  classification  is  shown  in  Figure  26. 


Figure  26.  Rule  image  for  trail  classified  pixels  generated  during  maximum  likelihood 
classification,  higher  values  represent  higher  likelihood  that  pixel  belongs  to 
trail  class.  The  upper  portion  of  the  image  shows  a  trail  segment,  and  the 
intermittent  long  diagonal  corresponds  to  a  stream  bed. 

Using  the  same  procedures,  several  classification  techniques  were  evaluated  and 
compared.  The  evaluation  and  comparison  will  be  further  explained  and  discussed  in 
Chapter  V.  ENVI  provides  nine  supervised  classification  techniques  that  have  been 
developed  and  used  primarily  for  spectral  classification  of  multispectral  and 
hyperspectral  imagery.  Of  these  nine,  four  were  chosen  and  evaluated  for  trail 
classification  on  the  Swanton  Pacific  Ranch  scene.  Maximum  likelihood  classification  as 
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discussed  in  Chapter  II  was  chosen  as  the  preferred  method.  Parallelepiped  classification 
uses  a  decision  rule  to  classify  data.  The  decision  boundaries  fonn  an  n-dimensional 
parallelepiped  classification  in  the  image  data  space  defined  based  on  the  standard 
deviation  threshold  from  the  mean  of  each  class  (ITT  Visual  Information  Solutions, 
2010).  If  the  pixels  in  the  scene  fall  within  the  threshold  it  belongs  to  that  class,  if  it  falls 
in  multiple  classes  it  is  assigned  to  the  last  class  matched,  and  if  it  does  not  fall  within  a 
threshold  it  is  unclassified.  The  minimum  distance  classifier  uses  mean  vectors  for  each 
training  set,  and  calculates  the  Euclidian  distance  from  each  unknown  pixel  to  the  mean 
vector  for  each  class  in  the  training  set.  All  unknown  pixels  are  then  classified  to  the 
closest  vector  from  the  training  set  classes.  The  Mahalanobis  distance  classification  uses 
statistics  for  each  class,  similar  to  maximum  likelihood  classification,  but  assumes  class 
covariances  are  the  same  so  it  is  a  faster  method. 

C.  POST  CLASSIFICATION  FILTER  METHODS 

Post  classification  image  processing  techniques  were  applied  to  classified  images 
to  improve  the  accuracy  of  the  final  product.  Convolution  and  morphology  filters  were 
investigated,  both  in  combination  and  individually.  ENVI  also  provides  sieve  and  clump 
techniques  that  can  be  applied  to  classified  image  products.  These  were  evaluated  as 
well.  The  output  results  from  each  of  the  filters  were  evaluated  and  will  be  discussed  in 
Chapter  V. 

The  convolution  filters  produce  images  in  which  the  brightness  value  of  a  given 
pixel  is  a  function  of  some  weighted  average  of  the  brightness  of  the  surrounding  pixels. 
The  extent  of  the  surrounding  pixels  considered  by  the  convolution  function  can  vary  in 
size,  and  is  known  as  a  kernel.  Median  and  Laplacian  convolution  filters  were  used  in 
this  research.  Median  filters  smooth  an  image,  removing  regions  of  noise  from  an  image 
smaller  than  the  size  of  a  user  specified  kernel.  The  Laplacian  filter  is  a  second 
derivative  edge  enhancement  filter  that  is  not  dependant  on  edge  direction. 

Morphological  operations  in  ENVI  are  dilation,  erosion,  opening,  and  closing. 
Dilate  fills  holes  smaller  than  the  user  selected  kernel  in  images,  where  erode  removes 
small  islands  of  pixels  that  are  smaller  than  the  kernel.  Opening  erodes  the  image 
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followed  by  dilation  by  the  same  size  kernel.  Opening  is  used  to  smooth  contours,  and 
remove  islands  and  peaks  in  an  image.  Closing  dilates  an  image  followed  by  erosion,  and 
is  used  to  fuse  narrow  breaks  and  fill  small  holes. 


^  1r  2  1  1  > 

(a)  (b)  (c)  (d)  (e)  (f) 


Figure  27.  Examples  of  binary  image  morphology:  (a)  original  image;  (b)  dilation;  (c) 
erosion;  (d)  majority;  (e)  opening;  (f)  closing.  Kernel  size  of  5x5  used  for 
each.  Majority  was  not  used  in  this  research,  but  rounds  sharp  comers 
(From  Szeliski,  2011) 
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V.  ANALYSIS 


A.  METHODS  OF  EVALUATION  AND  ANALYSIS 

To  evaluate  the  different  classification  techniques  and  filtering  methods  a  random 
sampling  of  pixels  from  the  scene  were  selected.  These  pixels  were  ground  truthed  and 
then  compared  to  classification  images  and  filtered  images  to  detennine  the  performance 
of  each  method.  This  evaluation  was  conducted  in  two  tiers.  The  first  tier  evaluated  the 
four  different  classification  techniques  used:  maximum  likelihood,  minimum  distance, 
parallelepiped,  and  Mahalanobis  distance.  After  evaluation,  maximum  likelihood  was 
chosen  as  the  best  technique,  and  post  classification  filtering  was  conducted  on  maximum 
likelihood  classifications  only.  The  second  tier  of  evaluations  compared  the  perfonnance 
of  different  filters  and  combinations  of  filters  to  each  other. 


To  detennine  the  number  of  pixels  needed  for  a  valid  sampling,  this  study  used 
the  margin  of  error  equation  (De  Veaux,  Velleman,  &  Bock,  2009): 

ME  =  z*  p 
V  n 


Where  z*  is  1.96  for  a  95%  confidence  level,  p=0.5  to  determine  the  largest  sample  size 
regardless  of  the  true  proportions,  and  ME  is  set  at  0.06.  The  equation  is  then: 


0.06  =  1.96, 


(0.5)(0.5) 


_  1.967(0.5X0.5) 
~  006 
n  =  (16. 3)2  =  266.8 


For  this  research  this  value  was  rounded  up,  and  300  pixels  were  randomly  selected  by 
the  ENVI  post  classification  random  sample  generator.  The  points  were  selected  from  a 
spatial  subset  of  Swanton  Pacific  Ranch  that  is  forested  and  has  many  roads  and  trails,  as 
seen  in  Figure  28. 
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Figure  28.  Spatial  subset  of  random  sampling.  Imagery  (on  left)  provided  by  NPS 
Remote  Sensing  Center  and  ENVI  slope  product  (on  right) 

Each  of  the  points  were  classified  as  trail  or  not  trail  using  a  combination  of 
ground  truth  survey  and  DEM  products.  The  sampling  followed  the  color  conventions 
from  the  classes  identified  in  the  training  set,  with  green  points  corresponding  to  trail 
classified  points  and  red,  blue,  and  yellow  corresponding  to  a  not  trail  classification.  Of 
the  267  points  generated,  after  ground  truth,  61  points  were  classified  as  trail  and  206 
were  classified  as  not  trail.  The  full  list  of  points  is  in  the  Appendix,  but  the  point 
breakdown  summary  was  as  follows: 


Table  3.  Random  point  sampling  summary 


ROI  name: 

ROI  rgb  value: 

ROI  npts: 

Random  Sample  (first_attempt_topol5July  /  [Green]  201  points) 

{0,  255,  0} 

117 

Random  Sample  (first_attempt_topol5July  /  [Red]  1228  points) 

{255,  0,  0} 

46 

Random  Sample  (first_attempt_topol5July  /  [Blue]  585  points) 

{0,  0,  255} 

65 

Random  Sample  (first_attempt_topol5July  /  [Yellow]  2616  points) 

{255,  255,  0} 
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Using  the  same  267  points,  each  of  the  classification  techniques  and  filtering 


methods  were  compared  to  the  ground  truth  to  evaluate  performance.  Perfonnance  was 
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measured  using  receiver  operator  characteristic  (ROC)  graphs  for  a  discrete  classifier. 
The  thematic  products  from  each  of  the  processes  in  this  research  result  in  a  classification 
into  one  of  only  two  classes.  With  a  discrete  two  class  classifier  and  an  instance  (known 
as  ground  truth),  there  are  four  possible  outcomes  (Fawcett,  2004).  If  the  instance  is 
positive  (trail)  and  it  is  classified  as  trail,  it  is  counted  as  a  true  positive;  if  it  is  negative 
(classified  not  trail),  it  is  counted  as  a  false  negative.  If  the  instance  is  negative  (not  trail) 
and  it  is  classified  as  negative  it  is  counted  as  a  true  negative;  if  it  is  classified  as  positive 
it  is  counted  as  a  false  positive.  Given  these  four  categories  a  2  x  2  confusion  matrix  can 
be  constructed  and  forms  the  basis  of  many  commonly  used  metrics  as  seen  in  Figure  29. 


Y 

Hypothesized 

class 

N 


Column  totals:  P  N 

fp  rate  —  tp  rate  =  ^ 

precision  =  tp+fp  recall  —  Lp- 

accuracy  =  F-measure  =  -r-, —  --1 

J  P+N  1  /  precision-f  1  /recall 

Figure  29.  Confusion  matrix  and  common  performance  metrics  calculated  from  it 
(From  Fawcett,  2004) 

Using  the  metrics  from  Figure  28,  ROC  graphs  are  plotted  with  the  TP  rate  on  the 
Y  axis  and  the  FP  rate  on  the  X  axis.  Generally,  one  point  in  ROC  space  is  better  than 
another  if  it  plots  in  the  upper  left  of  the  graph,  that  is,  it  has  a  high  true  positive  rate  and 
low  false  positive  rate.  Classifiers  on  the  left  side  of  the  graph  with  low  false  positive 


True  class 

P  n 


True 

Positives 

False 

Positives 

False 
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rates  make  positive  classifications  only  with  strong  evidence.  Alternatively  classifiers  on 
the  upper  right  side  of  the  graph  classify  nearly  all  positives  correctly,  but  have  high  false 
positive  rates.  A  diagonal  line  where  x=y,  corresponds  to  random  guessing. 

B.  EVALUATION  OF  CLASSIFICATION  TECHNIQUES 

Four  classification  techniques  were  evaluated  using  true  positive  rates  and  false 
positive  rates  as  measures  of  effectiveness.  Higher  true  positive  rates  and  lower  false 
positives  rates  are  better  for  each  discrete  classifier.  Each  of  the  four  classification 
techniques  rely  on  the  statistics  of  the  training  set  to  make  a  determination  as  to  which 
class  a  pixel  belongs,  as  discussed  in  Chapter  IV.  The  training  set  mean,  band  list  with 
description,  and  standard  deviations  are  provided  in  Tables  4  through  7. 


Table  4.  Training  set  mean  for  each  topographic  modeling  band  by  class 


0123456789  10  11 

Topographic  Modeling  Bands 
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Table  5.  Training  set  mean  for  topographic  modeling  bands  3  through  10  by  class 


ROI  Mean  (Bands  3  through  10) 


Table  6.  Topographic  Modeling  band  list 


Band 

Name 

Description 

1 

Slope 

Convention  ofO  degrees  for  a  horizontal  plane 

2 

Aspect 

The  direction  (azimuth)  that  a  surface  laces,  in  degrees  clockwise  from  North  (0  deg) 

3 

Shaded  Relief 

Renders  terrain  in  3D  by  use  of  shadows  that  would  be  cast  by  the  sun  from  the  NW 

4 

Profile  Convexity 

Rate  of  change  ofthe  slope  intersecting  with  the  plane  of  the  z-axis  and  aspect  direction 

5 

Plan  Convexity 

Rate  of  change  ofthe  aspect  intersecting  with  the  x,y  plane 

6 

Longitude  Convexity 

Measures  of  the  surface  curvature  orthogonally  in  the  down  slope  direction 

7 

Cross  Section  Convexity 

Measures  ofthe  surface  curvature  orthogonally  in  the  across  slope  direction 

8 

Minimum  C  urvature 

Minimum  overall  surface  curvature 

9 

Maximum  Curvature 

Maximum  overall  surface  curvature 

10 

RMS  Error 

Indication  of  how  well  the  quadratic  surface  fits  the  digital  elevation  model 

11 

Slope  Percentage 

The  percentage  or  degree  change  in  elevation  over  distance 

49 


Table  7.  Training  set  standard  deviation  for  all  topographic  modeling  bands  by  class 
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The  Mahalanobis  distance  classifier  had  the  lowest  false  positive  rate  of  the  four 


classifiers  with  23.7%.  It  ranked  three  out  of  four  for  true  positive  rate  with  67.2%. 
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Figure  30.  Mahalanobis  distance  confusion  matrix  with  metrics  and  thematic  map 
showing  trails  in  green 

The  maximum  likelihood  classifier  had  the  best  true  positive  rate  at  83.6%,  and 
ranked  two  out  of  four  for  false  positive  rate  with  3 1 .9%. 
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Figure  3 1 .  Maximum  likelihood  confusion  matrix  with  metrics  and  thematic  map 
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The  minimum  distance  classifier  had  the  second  highest  true  positive  rate  with 
83.6%,  and  also  had  the  highest  false  positive  rate  with  5 1.7%. 
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Figure  32.  Minimum  distance  confusion  matrix  with  metrics  and  thematic  map 


The  parallelepiped  classifier  had  the  lowest  true  positive  rate  with  63.9%,  and  had 
the  second  highest  false  positive  rate  with  37.2%. 
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Figure  33.  Parallelepiped  confusion  matrix  with  metrics  and  thematic  map 
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The  maximum  likelihood  classifier  was  chosen  as  the  best  classification  technique 
in  this  instance  due  to  its  high  true  positive  rate,  and  when  graphed  in  ROC  space  in 
Table  8.  The  Mahalanobis  distance  classifier  was  ranked  second  over  the  minimum 
distance  and  parallelepiped  classifiers  due  to  its  lower  false  positive  rate.  For  the 
remainder  of  the  research,  the  maximum  likelihood  classifier  results  were  used  to 
evaluate  different  filtering  methods. 


Table  8.  Graph  for  classification  techniques 


Graph  for  Classification  Techniques 
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C.  EVALUATION  OF  IMAGE  PROCESSING  METHODS 

The  same  evaluation  techniques  were  applied  to  the  resulting  products  from  the 
different  image  processing  methods.  Generally,  three  different  approaches  were  applied 
to  the  maximum  likelihood  classified  image.  The  first  approach  was  to  apply  a  median 
filter  to  remove  noise  in  the  image  followed  opening,  closing,  erosion,  and  dilation 
operators  of  varying  structural  size  and  comparing  the  results.  The  second  approach  was 
to  use  a  post  classification  tool  in  ENVI  named  sieve  and  clump.  The  sieve  function 
removes  isolated  pixels  using  blob  grouping  (ITT  Visual  Infonnation  Solutions,  2010). 
The  function  looks  at  local  neighborhood  for  each  pixel  and  determines  if  the  pixel  is 
grouped  with  pixels  of  the  same  class,  if  not  it  re-labels  the  pixel  unclassified.  The  clump 
function  groups  together  similarly  classified  areas  after  they  have  been  sieved.  The  third 
approach  used  a  decision  tree  to  create  a  final  product.  Using  the  decision  tree  the  first 
node  separates  trail  classified  points  from  non  trail  classified  points.  The  second  node  of 
the  decision  tree  takes  the  trail  classified  points,  known  as  node  one  survivors,  and 
applies  a  Laplacian  edge  detector.  The  output  from  this  approach  represents  the 
intersection  of  the  survivors  from  both  nodes  of  the  decision  tree  in  a  binary  mask. 

The  naming  conventions  used  attempt  to  capture  the  image  processing  method 
used  as  well  as  the  size  of  the  structuring  element.  This  was  done  in  an  effort  to  aid 
further  work,  for  example  an  image  that  has  a  median  filter  applied  with  a  5  x  5  kernel 
followed  by  the  morphological  operation  closing  with  a  3  x  3  kernel  would  be  labeled: 
median5by5_closing3by3.  The  sieve  and  clump  products  follow  a  similar  naming 
convention,  and  all  decision  tree  products  start  with  the  survivormask  and  follow  the 
naming  conventions  for  any  processing  that  was  applied  to  that  mask.  A  summary  of  the 
image  processing  techniques  used  and  their  respective  metrics  are  listed  in  Table  9. 


54 


Table  9.  Image  processing  techniques  used  with  perfonnance  metrics 


Method 

true  pos 

false  pos 

tru  neg 

false  neg 

tp  rate 

fp  rate 

precision 

accuracy 

maxjikelihood 

51 

66 

141 

9 

0.836 

0.319 

43.6% 

71.6% 

median5by5 

51 

66 

141 

9 

0.836 

0.319 

43.6% 

71.6% 

median5by5 opening3by3 

51 

66 

141 

9 

0.836 

0.319 

43.6% 

71.6% 

median5by5 closing3by3 

50 

66 

141 

10 

0.820 

0.319 

43.1% 

71.3% 

median5by5 closing5by5 

41 

65 

142 

19 

0.672 

0.314 

38.7% 

68.3% 

median5by5 erode3by3 dilate7by7 

32 

64 

143 

28 

0.525 

0.309 

33.3% 

65.3% 

median5by5 erode7by7 dilate3by3 

48 

71 

136 

12 

0.787 

0.343 

40.3% 

68.7% 

sieveAllclasses clumpAllclasses 

52 

66 

141 

8 

0.852 

0.319 

44.1% 

72.0% 

sieveAllclasses clumpgeen 

52 

67 

140 

8 

0.852 

0.324 

43.7% 

71.6% 

survivormask 

15 

9 

198 

45 

0.246 

0.043 

62.5% 

79.5% 

survivormask median3by3 

14 

4 

226 

56 

0.230 

0.019 

77.8% 

89.6% 

survivormask_median3by3_3cyclesclosing 

20 

7 

223 

50 

0.328 

0.034 

74.1% 

90.7% 

The  image  processing  techniques  only  increased  true  positive  rates  in  the  sieve 
then  clump  processes,  however,  only  slightly  correctly  identifying  only  one  additional 
trail  point  than  the  original  classifier.  The  most  accurate  processes  were  the  products 
from  the  decision  tree.  All  three  decision  tree  products  reduced  the  false  positive  rate  to 
less  than  5%,  with  the  best  performer  with  regards  to  accuracy  being  the  decision  tree 
product,  followed  by  a  median  3x3  filter,  followed  by  3  cycles  of  closing  with  a  3  x  3 
filter.  This  process  works  the  best  since  it  moves  from  the  more  liberal  classification 
process  and  step-by-step  removes  unwanted  artifacts  based  on  trail  characteristics.  The 
maximum  likelihood  classification  identifies  a  high  percentage  of  trails  in  the  scene  with 
an  83.6%  true  positive  rate.  The  Laplacian  convolution  enhances  edges  of  long  linear 
features  (trails)  without  regard  to  direction,  and  the  decision  tree  survivors  are  those 
pixels  that  are  both  classified  as  trail  and  have  edges.  This  product  is  smoothed  with  a 
median  filter  to  remove  noise  in  the  scene  while  preserving  edges,  and  finally  the  closing 
fuses  any  narrow  breaks  and  fills  small  holes  in  the  trail  segments.  The  thematic  map 
from  this  process  is  shown  in  Figure  34. 
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Figure  34.  Binary  image  showing  trail  classified  points  in  black  from  image  processing 


Figure  35. 


VI.  SUMMARY 


The  overall  result  of  this  experiment  represented  by  the  analysis  in  Chapter  V  is 
very  encouraging  that  the  process  of  statistical  classification  followed  by  image 
processing  can  correctly  identify  roads  and  trails  under  canopy  using  LiDAR.  The  point 
densities  of  current  LiDAR  systems  are  capable  of  penetrating  second  growth  forest 
canopy,  and  producing  accurate  digital  terrain  models  (DTMs).  In  this  data  set  the  slope 
values,  convexities  along  different  planes,  and  curvatures  provided  sufficient  statistics  to 
characterize  unpaved  road  and  trail  segments  within  the  study  area.  Supervised 
classification  techniques  traditionally  used  in  remote  sensing  on  multispectral  and 
hyperspectral  images  were  applied  successfully  to  identify  roads  and  trails  after  training 
data  was  identified. 

Four  classification  techniques  were  used  in  the  experiment  to  identify  the  best 
techniques  for  classifying  roads  and  trails  in  forested  areas.  Among  maximum 
likelihood,  Mahalanobis  distance,  minimum  distance,  and  parallelepiped  classifiers 
maximum  likelihood  produced  the  highest  true  positive  rate  and  was  chosen  as  the  sole 
classifier  to  conduct  image  processing.  However,  the  Mahalanobis  distance  classifier  did 
result  in  the  highest  accuracy  rating  and  may  warrant  more  investigation  in  further  work. 

Image  processing  techniques  were  applied  to  the  result  of  the  maximum 
likelihood  classification  and  compared  to  each  other.  Three  different  approaches  were 
used  and  evaluated  against  each  other.  The  sieve  and  clump  operators  maintained  the 
highest  true  positive  rate.  The  results  from  the  decision  tree  approach  reduced  the  true 
positive  rate,  but  also  had  the  lowest  false  positive  rate  which  led  to  it  having  the  highest 
accuracy  rate. 

In  the  evaluations  of  both  the  classification  techniques  and  the  image  processing 
methods,  the  best  techniques  depend  on  the  application.  The  maximum  likelihood 
classification  and  the  sieve  and  clump  processing  can  be  thought  of  as  more  “liberal”: 
they  classify  nearly  all  positives  correctly,  but  have  a  high  false  positive  rate.  The 
Mahalanobis  distance  classifier  and  the  decision  tree  processing  techniques  are  more 
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“conservative”:  they  make  few  false  positive  classifications,  but  their  true  positive  rate  is 
low  as  well.  Either  set  of  techniques  results  in  a  thematic  map  that  successfully  maps 
roads  and  trails  under  canopy. 
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VII.  CONCLUSION 


The  results  of  the  experiment  conducted  for  this  research  indicate  that  roads  and 
trails  can  be  automatically  identified  under  dense  tree  canopy  using  LiDAR  data.  Using 
classification  techniques  alone  true  positive  rates  of  over  83%  can  be  achieved, 
identifying  trails  under  dense  redwood  and  Douglas  fir  second  growth  forests.  Applying 
image  processing  techniques  to  this  classification  product  can  reduce  false  positive  rates 
to  below  4%,  increasing  the  accuracy  of  the  trail  identification  to  90.7%.  Thematic  map 
products  from  these  processes  could  prove  invaluable  for  both  military  and  commercial 
applications.  During  intelligence  preparation  of  the  battlefield  (IPB),  commanders  and 
intelligence  professionals  can  use  this  process  to  define  the  lines  of  communications 
(LOC),  which  an  enemy  threat  may  use  to  move  personnel  and  supplies  through  an  area 
of  operations. 

In  the  field  of  LiDAR,  there  are  many  ongoing  research  and  development 
programs  in  government  and  commercial  sponsored  studies.  Areas  of  follow-on  research 
related  to  this  research  include: 

•  Further  characterize  trails  to  increase  true  positive  rates  and  reduce  false 
positive  rates. 

•  Test  this  process  against  other  data  sets  with  varying  terrain  and  point 
densities. 

•  Develop  geometric  pathway  extraction  algorithm  to  produce  centerline 
and  trail  width  vectors. 

•  Write  code  to  fully  automate  steps  in  this  process  from  point  cloud  to 
refined  trail  map  output. 

The  accuracy  of  road  and  trail  network  mapping  using  LiDAR  can  be  used  for 
quantitative  terrain  analysis  without  the  need  for  ground  reconnaissance  in  areas 
unobservable  to  electro-optical  imagery.  LiDAR  provides  the  ability  to  determine  road 
networks  on  large  scales  in  denied  areas  where  ground  survey  is  not  available.  This 
research  demonstrates  the  value  of  LiDAR  collected  data  in  areas  where  traditional 
remote  sensing  techniques  for  intelligence  preparation  of  the  battlefield  are  insufficient. 
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LiDAR  holds  great  potential  to  provide  accurate,  detailed,  and  large  coverage  area 
support  to  ground  maneuver  forces  deployed  in  diverse  and  complex  environments. 
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APPENDIX  -  RANDOM  POINT  SAMPLING  LISTS 


ROI  name:  Random  Sample  (first_attempt_topol5July  /[Green]  201  points) 


ROI  rgb  value;  {0,  255,  0} 


ROI  ripts:  117 


Map  X 


606310S.62 


6061898.27 


6063853.65 


6062994.07 


6063604.3 


6062803.78 


6061845.78 


6065011.78 


6062725.04 


6064864.14 


6063131.86 


6062803.78 


6063079.37 


6062692.23 


6061727.67 


6061635.8 


Map  Y 


1850271.04 


1850251.35 


1850248.07 


1850244.79 


1850244.79 


1850205.42 


1850195.58 


1850189.02 


1850179.17 


1850175.89 


1850162.77 


1850152.93 


1850133.24 


1850080.75 


1850047.94 


1850041.38 


Lat 


37.06384 


37.06373 


37.06382 


37.06377 


37.0638 


37.06365 


37.063S7 


37.06371 


37.06357 


37.06367 


37.063S5 


37.0635 


37.06346 


37  0633 


37.06316 


37.06314 


6062692.23 

1850024.98 

37.06315 

6062902.2 

1850021.69 

37.06315 

6065320.18 

1849992.17 

37.06319 

6062505.22 

1849979.04 

37.06301 

6061986.85 


6061356.93 


6063469.79 


6061484.88 


6063394.33 


6063660.08 


6063174.51 


6063187.64 


606S149.S7 


6061409.42 


6061501.29 


6061304.44 


6061419.27 


6061478.32 


6061386.46 


6061963.89 


6062974.38 


6064312.96 


6061120.71 


606214105 


6061937.64 


6062269 


6062104.96 


6061534.1 


6061596.43 


6063276.22 


1849962.64 


1849939.67 


1849870.78 


1849854.37 


1849811.72 


1849769.07 


1849749.39 


1849746.1 


1849723.14 


1849716.58 


1849713.3 


1849680.49 


1849673.93 


1849667.36 


1849611.59 


1849611.59 


1849598.47 


1849575.5 


1849532.85 


1849434.43 


1849431.14 


1849427.86 


1849427.86 


1849427.86 


1849395.06 


1849316.32 


37.06294 


37.06285 


37.06276 


37.06262 


37.0626 


37.06249 


37.06241 


37.06241 


37.06244 


37.06224 


37.06223 


37.06213 


37.06212 


37.0621 


37.06195 


37.06198 


37.06199 


37.06199 


37.06172 


37.0615 


37.06148 


37.06149 


37.06148 


37.06145 


37.06136 


37.06123 


Lon 


-122.209 


-122.213 


-122.206 


122.209 


-122.207 


-122.21 


122.213 


122.202 


-122.21 


122.203 


-122.209 


-122.21 


122.209 


-122  21 


-122.213 


122.214 


-122.21 


-122.209 


122.201 


122.211 


122.213 


-122.215 


122.207 


-122.214 


122.208 


-122.207 


122.208 


-122.208 


-122.202 


122.215 


122.214 


122.215 


122.214 


-122.214 


122.215 


122.213 


122.209 


122.205 


122.215 


122.212 


-122.213 


122.212 


-122.212 


122.214 


-122.214 


122.208 


61 


62 


98 

1536 

2795 

6062928.45 

1847406.87 

37.05597 

-122.209 

1 

99 

1601 

2809 

6063141.7 

1847360.94 

37.05585 

-122.208 

1 

100 

1221 

2809 

6061894  99 

1847360  94 

37  05579 

-122.213 

1 

101 

2096 

2814 

6064765.72 

1847344.53 

37.05589 

-122,203 

1 

102 

997 

2817 

6061160.08 

1847334.69 

37.05568 

-122.215 

1 

103 

2064 

2817 

6064660.73 

1847334.69 

37.05586 

-122.203 

1 

104 

2145 

2818 

6064926.48 

1847331.41 

37.05586 

-122.202 

1 

105 

1707 

2820 

6063489.47 

1847324.85 

37.05577 

-122.207 

1 

106 

1031 

2831 

6061271.63 

1847288.76 

37.05556 

-122.215 

1 

107 

1436 

2842 

6062600.37 

1847252.67 

37.05553 

122.21 

1 

108 

1178 

2844 

6061753.91 

1847246.11 

37.05547 

-122.213 

1 

109 

2179 

2851 

6065038.03 

1847223.14 

37.05557 

-122.202 

1 

110 

1684 

2855 

6063414.01 

1847210.02 

37.05545 

-122.207 

1 

111 

1752 

2857 

6063637.11 

1847203.46 

37.05545 

-122.207 

1 

112 

1022 

2857 

6061242.1 

1847203.46 

37.05533 

-122.215 

1 

113 

1260 

2870 

6062022.94 

1847160.81 

37.05525 

-122.212 

1 

114 

2254 

2870 

6065284.09 

1847160.81 

37.05541 

-122.201 

1 

115 

1984 

2872 

6064398.26 

1847154.25 

37  05535 

-122  204 

1 

116 

2107 

2875 

6064801.81 

1847144.4 

37.05534 

-122.203 

1 

117 

1638 

2879 

6063263.1 

1847131.28 

37.05523 

-122.208 

1 

63 
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ROI  name:  Random  Sample  (fir$t  auempt_topol5Julv  /I Blue]  5S5  points) 


ROIrgb  value:  (0,0,  255) 


ROI  npts:  65 


Map  X 


1721  6063335.3 


1748  6065349.7 


1757  6062689 


1825  6063292.6 


1843  6065503.9 


1875  6062692.2 


1885  6062639.7 


1891  6063696.2 


1910  6061704.7 


1943  6062436.3 


2034  6065579.4 


2052  6065274.3 


2076  6064677.1 


2092  6064319.5 


2096  6064621.4 


2105  6064965.9 


2131  6061721.1 


2133  6065103.6 


2141  6061678.5 


2165  6065244.7 


2226  6064296.6 


2227  6062452.7 


2293  6064306.4 


2301  6061708 


2306  6061714.5 


2328  6061635.8 


2332  6064641.1 


2341  6064421.2 


2344  6064411.4 


2403  6062367.4 


2411  6065530.2 


2418  6061534.1 


2442  6061432.4 


2489  6061064.9 


2494  6062790.7 


2569  6064132.5 


2582  6062879.2 


2592  6065175.8 


2598  6062879.2 


2600  6062754.6 


2608  6062649.6 


2626  6064076.7 


2626  6062918.6 


2634  6062741.4 


2646  6062967.8 


2652  6064988.8 


2660  6065287.4 


2669  6064358.9 


Map  Y 


1850930.5 


1850841.9 


1850812.4 


1850589.3 


1850530.2 


1850425.2 


1850392.4 


1850372.7 


1850310.4 


18S0202.1 


1849903.6 


1849844.5 


1849765.8 


1849713.3 


1849700.2 


1849670.7 


1849585.3 


1849578.8 


1849473.8 


1849273.7 


1849270.4 


1849053  9 


1849027.6 


1849011.2 


1848939 


1848925.9 


1848896.4 


1848693 


1848666.7 


1848643.7 


1848565 


1848410.8 


18483944 


1848148.3 


1848105.7 


1848072.9 


1848053.2 


1848046.6 


1848020.4 


1847961.3 


1847961.3 


1847935.1 


1847895.7 


1847876 


1847849.8 


1847820.3 


Lat 


37  06567 


37.06552 


37.06531 


37.06473 


37.06468 


37.06425 


37.06415 


37.06415 


37.06388 


37.06362 


37  06296 


37.06278 


37.06254 


37.06237 


37.06235 


37.06229 


37.06189 


37.06204 


37.0618 


37.06176 


37.06116 


37  06106 


37  06056 


37.06036 


37  06031 


37.06011 


37  06023 


37  06013 


37.06011 


37.05947 


37.05956 


37.0593 


37.05908 


37.05863 


37.05867 


37  05807 


37.05789 


37  05791 


37.05774 


37  05772 


37.05764 


37.05755 


37.05749 


37.05741 


37.05731 


37.05736 


37.0573 


37.05718 


Lon 


-122.208 


122.2011 


122.2102 


-122.2081 


-122.2005 


122.2102 


122.2103 


122.2067 


122.2135 


-122.211 


122.2002 


122.2013 


-122.2033 


-122.2045 


-122.2035 


-122.2023 


-122.2134 


-122.2018 


122.2136 


-122.2014 


-122.2046 


-122.2109 


-122.2045 


-122.2134 


122.2134 


-122.2137 


-122.2034 


-122.2041 


122.2042 


-122.2112 


122.2003 


-122.214 


122.2144 


-122.2156 


122.2097 


-122.2051 


-122.2094 


122.2015 


-122.2094 


122.2098 


122.2102 


122.2053 


-122.2092 


-122.2098 


-122.2091 


122.2021 


-122.2011 


-122.2043 


Not  Trail  =  1 


65 


49 

1275 

2674 

6062072.2 

1847803.9 

37.05702 

-122.2121 

1 

50 

1123 

2686 

6061573.5 

1847764.5 

37.05688 

122.2138 

1 

51 

2303 

2690 

6065444.9 

1847751.4 

37.05704 

122.2006 

1 

52 

1641 

2699 

6063272.9 

1847721.8 

37.05685 

-122.208 

1 

53 

1494 

2702 

6062790.7 

1847712 

37.0568 

-122.2097 

1 

54 

1897 

2716 

6064112.8 

1847666.1 

37.05674 

-122.2051 

1 

55 

1686 

2725 

6063420.6 

1847636.5 

37.05663 

-122.2075 

1 

56 

1171 

2726 

6061731 

1847633.3 

37.05653 

-122.2133 

1 

57 

2285 

2750 

6065385.8 

1847554.5 

37.0565 

-122.2008 

1 

58 

1021 

2769 

6061238.8 

1847492.2 

37.05612 

-122.215 

1 

59 

1855 

2783 

6063975 

1847446.2 

37.05613 

-122.2056 

1 

60 

1466 

2808 

6062698.8 

1847364.2 

37.05584 

-122.2099 

1 

61 

1287 

2811 

6062111.5 

1847354.4 

37.05578 

-122.212 

1 

62 

1199 

2838 

6061822.8 

1847265.8 

37.05553 

-122.2129 

1 

63 

1093 

2840 

6061475 

1847259.2 

37.05549 

122.2141 

1 

64 

1119 

28S2 

6061560.3 

18472195 

37.05539 

-122.2138 

1 

65 

1080 

2892 

6061432.4 

1847088.6 

37.05502 

-122.2143 

1 

66 
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